############################################################
#### HOW UNCERTAINTY ABOUT WAR OUTCOMES AFFECTS WAR ONSET 
#### BAS AND SCHUB
#### JCR REPLICATION FILE 
############################################################


######################
## CONTENTS OVERVIEW
######################

# (I) Figure 1 from paper

# (II) Table 1 from paper

# (III) Table 2 from Online Appendix

# (IV) Table 7 from Online Appendix


########
# Set working directory and load packages
########
library(foreign)
library(lme4)


#########################
# (I) Figure 1 from paper
#########################
fig<-read.dta("SystemYear.dta")

plot(fig$year, fig$W_SystemRELa,type="l",lwd=4, ylab="System Uncertainty",xlab="Year",main="System Uncertainty Over Time")
points(fig $year, fig$SystemRELa,type="l",col="grey55",lwd=4)
legend("topright",legend=c("Weighted","Unweighted"),col=c("black","grey55"),lwd=c(4,4),bty="n")

#########################
# (II) Table 1 from paper
#########################

data<-read.dta("DyadYear.dta")
head(data)

########
## (A) Lower portion ##

#i) Top 6 state cinc scores 1890
rev(sort(unique(c(data$cap_1[data$year==1890],data$cap_2[data$year==1890]))))[1:6]

#ii) Top 6 state cinc scores 1992
rev(sort(unique(c(data$cap_1[data$year==1992],data$cap_2[data$year==1992]))))[1:6]

#iii) Uncertainty and Dispersion scores 1890
data$SystemRELa[data$year==1890][1]
data$con[data$year==1890][1]

#iv) Uncertainty and Dispersion scores 1992
data$SystemRELa[data$year==1992][1]
data$con[data$year==1992][1]

#########
## (B) Upper portion ##

# creating the matrices

# 1855
casea <- matrix(c(0,.36, .32, .21, .18, .36, 0, .46, .32, .28, .32, .46, 0, .36, .31, .21, .32, .36, 0, .45, .18, .28, .31, .45, 0), 5, 5)
# making sure they are symmetric
casea-t(casea)
# 1912
caseb <- matrix(c(0,.38, .36, .34, .34, .38, 0, .47, .46, .45, .36, .47, 0, .48, .48, .34, .46, .48, 0, .49, .34, .45, .48, .49, 0), 5, 5)
caseb-t(caseb)
# 1955
casec <- matrix(c(0,.40, .26, .16, .14, .40, 0, .34, .22, .20, .26, .34, 0, .36, .33, .16, .22, .36, 0, .47, .14, .20, .33, .47, 0), 5, 5)
casec-t(casec)
# 1990
cased <- matrix(c(0,.47, .43, .30, .29, .47, 0, .46, .32, .31, .43, .46, 0, .36, .34, .30, .32, .36, 0, .49, .29, .31, .34, .49, 0), 5, 5)
cased-t(cased)

tmin <- matrix(0, 5,5)
tmax <- matrix(.5, 5,5)

# putting zeros in the diagonal
diag(tmax) <- 0

etmin <- eigen(tmin)$values[1]
etmax <- eigen(tmax)$values[1]
ecasea <- (eigen(casea)$values[1] - etmin)/(etmax-etmin)
ecaseb <- (eigen(caseb)$values[1] - etmin)/(etmax-etmin)
ecasec <- (eigen(casec)$values[1] - etmin)/(etmax-etmin)
ecased <- (eigen(cased)$values[1] - etmin)/(etmax-etmin)

# Uncertainty Scores
print(c(ecasea, ecaseb, ecasec, ecased))


#########################
# (III) Table 2 from Online Appendix
#########################

########
## (A) Load data
data<-read.dta("DyadYear.dta")


########
## (B) Produce Models 1 and 2 from Table 2 of Online Appendix

##### System Uncertianty Unweighted  #####
base<-lmer(neworigforce ~ SystemRELa + relcap + rival + contig + jdem + alliance + log(pceyrs+1) + log(pceyrs2+1) + log(pceyrs3+1) + (1|dyad) + (1|year), data=data, family=binomial)
summary(base)


##### System Uncertianty Weighted (by just capabilities)  #####
weight<-lmer(neworigforce ~ W_SystemRELa + relcap + rival + contig + jdem + alliance + log(pceyrs+1) + log(pceyrs2+1) + log(pceyrs3+1) + (1|dyad) + (1|year), data=data, family=binomial)
summary(weight)



#########################
# (IV) Table 7 from Online Appendix
#########################

##########
## (A) Load data
rjs<-read.dta("DyadYear.dta")
head(rjs)


##########
## (B) Run each model to be compared

#i) System Uncertainty (Model 1 of Table 7 in Online Appendix)
system<-glm(rjs$neworigforce ~ rjs$W_SystemRELa + rjs$relcap + rjs$rival + rjs$contig + rjs$jdem + rjs$alliance + rjs$pceyrs + rjs$pceyrs2 + rjs$pceyrs3, family=binomial (link=logit))
summary(system)

#ii) Dispersion (Model 2 of Table 7 in Online Appendix)
sysscore<-glm(rjs$neworigforce ~ rjs$con + rjs$relcap + rjs$rival + rjs$contig + rjs$jdem + rjs$alliance + rjs$pceyrs + rjs$pceyrs2 + rjs$pceyrs3, family=binomial (link=logit))
summary(sysscore)

#iii) Polarity (Model 3 of Table 7 in Online Appendix)
rjs$multi<-ifelse(rjs$year<1946,1,0)
rjs$bi<-ifelse(rjs$year>1945 & rjs$year<1992,1,0)

polarity<-glm(rjs$neworigforce ~ rjs$multi +rjs$bi + rjs$relcap + rjs$rival + rjs$contig + rjs$jdem + rjs$alliance + rjs$pceyrs + rjs$pceyrs2 + rjs$pceyrs3, family=binomial (link=logit))
summary(polarity)

#iv) Regional Uncertainty (Model 4 of Table 7 in Online Appendix)
region<-glm(rjs$neworigforce ~ rjs$W_Regionalb + rjs$relcap + rjs$rival + rjs$contig + rjs$jdem + rjs$alliance + rjs$pceyrs + rjs$pceyrs2 + rjs$pceyrs3, family=binomial (link=logit))
summary(region)


##########
## (C) Run Vuong test (alter z to change comparison model)

x<-system
z<-sysscore  #NOTE: alter between sysscore, polarity, and region
    y <- x$y
    N <- length(y)
    fam.1 <- family(x)$family
    fam.2 <- family(z)$family                                                              
    fit.mod.1 <- fitted(x)
    fit.mod.2 <- fitted(z)
    xb.1 <- predict.glm(x)
    xb.2 <- predict.glm(z)
    loglik.1 <- switch(fam.1,
        poisson = (xb.1*y-exp(xb.1)-log(gamma(y+1))),
        binomial = y*log(fit.mod.1)+(1-y)*log(1-fit.mod.1),
        gaussian = (-log(2*pi*sum((residuals(x))^2)/N))/2 - 
            (1/2)*((residuals(x))/sqrt(sum((residuals(x))^2)/N))^2)
    loglik.2 <- switch(fam.2,
        poisson = (xb.2*y-exp(xb.2)-log(gamma(y+1))),
        binomial = y*log(fit.mod.2)+(1-y)*log(1-fit.mod.2),
        gaussian = (-log(2*pi*sum((residuals(z))^2)/N))/2 - 
            (1/2)*((residuals(z))/sqrt(sum((residuals(z))^2)/N))^2)
    pdim <- x$rank
    if (fam.1 %in% c("gaussian", "Gamma", "inverse.gaussian"))
        pdim <- pdim+1
    qdim <- z$rank
    if (fam.2 %in% c("gaussian", "Gamma", "inverse.gaussian"))
        qdim <- qdim+1
    loglik.1.sum <- pdim-AIC(x)/2
    loglik.2.sum <- qdim-AIC(z)/2
    lldiff <- loglik.1-loglik.2
    vuong.den <- (sqrt(mean(lldiff^2)-((mean(lldiff))^2)))*sqrt(N)
    vuong.num <- (loglik.1.sum-loglik.2.sum-
        (((pdim/2)*log(N))-((qdim/2)*log(N))))
    Vuong <- vuong.num/vuong.den
    Vuong.p <- 2*(1-pnorm(abs(Vuong)))
c(Vuong,Vuong.p)

# System Uncertainty vs. Polarity: p<0.001
# System Uncertainty vs. Dispersion: p<0.01
# System Uncertainty vs. Regional Uncertainty: p<0.01






